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ABSTRACT 

Context. Twenty- six high-luminosity IRAS sources believed to be collection of stars in the early phases of high-mass star formation 
have been observed in the Near-IR (J, H, Kg) to characterize the clustering properties of their young stellar population and compare 
them with those of more evolved objects (e.g., Herbig Ae/Be stars) of comparable mass. All the observed sources possess strong con- 
tinuum and/or line emission in the millimeter, being therefore associated with gas and dust envelopes. Nine sources have far-IR colors 
characteristic of UCHII regions while the other 17 are likely being experiencing an evolutionary phase that precedes the Hot-Cores, 
as suggested by a variety of evidence collected in the past decade. 

Aims. To gain insight into the initial conditions of star formation in these clusters (Initial Mass Function [IMF], Star Formation History 
[SFH]), and to deduce mean values for cluster ages. 

Methods. For each cluster we carry out aperture photometry. We derive stellar density profiles, color-color and color-magnitude dia- 
grams, and color (HKCF) and luminosity (KLF) functions. These two functions are compared with simulated KLFs and HKCFs from 
a model that generates populations of synthetic clusters starting from assumptions on the IMF, the SFH, and the Pre-MS evolution, 
and using the average properties of the observed clusters as boundary conditions (bolometric luminosity, dust distribution, infrared 
excess, extinction). 

Results. Twenty-two sources show evidence of clustering with a stellar richness indicator that varies from a few up to several tens 
of objects, and a median cluster radius of 0.7 pc. A considerable number of cluster members present an infrared excess characteristic 
of young Pre-Main-Sequence objects. For a subset of 9 detected clusters, we could perform a statistically significant comparison of 
the observed KLFs with those resulting from synthetic cluster models; for these clusters we find that the median stellar age ranges 
between 2.5 • 10^ and 5-10^ years, with evidence of an age spread of the same entity within each cluster. We also find evidence that 
older clusters tend to be smaller in size, in line with the fact that our clusters are on average larger than those around relatively older 
Herbig Ae/Be stars. Our models allow us to explore the relationship of the mass of the most massive star in the cluster with both the 
clusters richness and their total stellar mass. Although such relationships are predicted by several classes of cluster formation models, 
their detailed analysis suggests that our modeled clusters may not be consistent with them resulting from random sampling of the 
IMF 

Conclusions. Our results are consistent with a star formation which takes place continuously over a period of time which is longer 
than a typical crossing time. 

Key words. Stars: formation - Stars: imaging - Stars: luminosity function, mass function - Stars: pre-main sequence - Infrared: stars 



1. Introduction 

The last few decades have been characterized by a large effort to 
improve our understanding of how stars form both from a theo- 
retical and from an observational point of view. As a result, today 
we have reached a good understanding of how isolated low-mass 
stars form (Klein at al. 2006). The widely accepted scenario is 
that low mass stars form through the gravitational collapse of a 
prestellar core followed at later stages by disk accretion. 

Extending this theory to high-mass stars is not trivial. High 
mass (proto-)stars reach the Zero Age Main Sequence while still 
accreting. When the central protostar reaches about 10 M© hy- 
drogen fusion ignites in the core and the star's radiation pressure 
and wind should prevent further accretion. This obviously is a 
paradox given that more massive stars do form. Several theo- 
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ries have been put forward to solve this dilemma (Zinnecker & 
Yorke l20()7]) : accretion rates as high as three orders of magnitude 
larger than in the case of low-mass stars (Cesaroni I2005K and 
non- spherical accretion geometries (Nakano'1989! Yorke '2002) 
Keto I2003K or coalescence in dense (proto-)stellar clusters 
(Bonnell etal. 179981 . 

All these theories have predictions that, in principle, could 
be tested observationally. In the last decade a large effort was 
made in trying to detect massive accretion disks (Cesaroni et 
al. i2006j , powerful outflows (Beuther et al. .2002, Cesaroni et 
al. 12005b and dense protostellar clusters (Testi et al. 119991 de 
Wit et al.'2005 ), all of these phenomena are predicted by one or 
the other formation theory. None of these efforts have provided 
conclusive arguments in favour or against any of the theories. 

In this paper we explore the properties of embedded clusters 
associated with high-mass protostellar candidates. Our sample 
was selected from a larger sample of candidate high-mass pro- 
tostars selected and analyzed, in the past decade, by Molinari et 
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al. ([T996l[T998l [20001 [20021) and Brand et al. ([2QQT1) . In Sect.|2] 

we present the observations and data analysis (source extraction, 
photometry), in Sect. [3] we discuss data elaboration and interpre- 
tation. In Sect. |4] we present our Synthetic Cluster Generation 
model and the method of comparison between synthetic and ob- 
served clusters and the results of using this technique. Finally 
in Sect. [5] we compare our objects with more evolved ones and 
present our conclusions. 

2. Observations and Data Analysis 

Program fields are listed in Table [T] and were imaged in J, H, and 
Ks bands. A total of 15 fields were observed in three nights in 
November 1998 at the Palomar 60-inch telescope equipped with 
a 256x256 NICMOS-3 array, with a pixel scale of 0:'62/pix and 
total FOV of 2f6x2.'6. The remaining 1 1 fields were observed in 
3 nights in August 2000 at the ESO-NTT using the 1024x1024 
SOFI camera, with a pixel scale of Of'29/pix and a total FOV 
of 4.'9x4f9. Standard dithering techniques were used to mini- 
mize impact of bad pixels and optimize flat-fielding, allowing to 
achieve for each field a total of 5min integration time per band 
(in the central portion of the observed field) for a covered area 
of 3f5x3.'5 for Palomar observations, and 20min (lOmin for the 
band Kg ) at NTT with a total covered area of 6f5x6.'5. Suitable 
calibration sources from the list of Hunt et al. ( 1988 ) were ob- 
served regularly during the observations to track atmospheric 
variations for diff'erent airmasses. Standard stars and target fields 
were observed at airmasses no greater than 1.7 at NTT, and 1.3 
at Palomar; we determined average zero-point magnitudes for 
each night and used them to calibrate our photometry. For each 
field the images in the three bands were registered and astromet- 
ric solutions were determined using a few bright optically visible 
sources. 

The images for all observed fields, with superimposed 
submillimeter continuum emission distribution when available 
(Molinari et al. 2008a) are pr esented in A ppendix [A| and avail- 
able online at http://galatea.ifsi-roma.inaf.it/faustini/K-images 

2.1. Point Source Extraction and Photometry 

The extraction and photometry of point sources for all images 
was carried using the IRAF package. The r.m.s of the back- 
ground signal and the FWHM of point sources were measured 
throughout the images to characterize the image noise and PSF 
properties; these parameters were fed to the DAOFIND task for 
source extraction, where a detection threshold of 3cr was used 
for all images. Sources with saturated pixels were excluded from 
the analysis; the linearity of the system response was checked a 
posteriori comparing, both for the Palomar and the NTT data, 
the magnitudes obtained to those from 2MASS using a few stars 
with magnitudes reaching up to the maximum values found in 
our photometry files; the relations between the 2MASS magni- 
tudes and ours in the three bands were found to be linear over 
the entire magnitude range of the detected sources. There were 
clearly brighter objects in the various fields, but their peaks were 
already flagged as saturated and were excluded from the detec- 
tion process. 

The photometry of sources is critical in very dense stellar 
fields like the inner Galactic Plane, where all of our target fields 
lie, and where the crowding is such that more than one source 
can enter any plausible aperture that can be chosen, or any an- 
nulus used for background estimation. This problem is of course 
magnified in the clustered environments that are detected in sites 
of massive star formation (see |3.1| below). 



Table 1. 
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Source 


IRAS Name 


Qf(J2000) 


(5(J2000) 


Tel. 


Mol^ 










3 


00420+5530 


00:44:57.6 


+55:46:52 


Pal 


8 


05137+3919 


05:17:13.3 


+39:22:14 


Pal 


9 


05168+3634 


05:20:16.2 


+36:37:21 


Pal 


11 


05345+3157 


05:37:47.8 


+31:59:24 


Pal 


12 


05373+2349 


05:40:24.4 


+23:50:54 


NTT 


15 


06056+2131 


06:08:41.0 


+21:31:01 


Pal 


28 


06584-0852 


07:00:51.0 


-08:56:29 


Pal 


30 


17450-1742 


17:48:09.3 


-27:43:21 


NTT 


38 


18024-2119 


18:06:18.0 


-21:42:00 


NTT 


45 


18144-1723 


18:17:24.2 


-17:22:13 


NTT 


50 


18162-1612 


18:19:07.5 


-16:11:21 


NTT 


59 


18278-1009 


18:30:35.2 


-10:07:12 


Pal 


75 


18511+0146 


18:53:38.1 


+01:50:27 


Pal 


82 


18565+0349 


18:59:03.4 


+03:53:22 


NTT 


84 


18567+0700 


18:59:13.6 


+07:04:47 


NTT 


98 


19092+0841 


19:11:37.4 


+08:46:30.0 


NTT 


99 


19094+0944 


19:11:52.0 


+09:49:46 


Pal 


103 


19213+1723 


19:23:37.0 


+ 17:28:59 


NTT 


109 


19374+2352 


19:39:33.2 


+23:59:55 


NTT^ 


110 


19388+2357 


19:40:59.4 


+24:04:39 


NTT^ 


136 


21307+5049 


21:32:31.5 


+51:02:22.0 


Pal 


139 


21519+5613 


21:53:38.8 


+56:27:53.0 


Pal 


143 


22172+5549 


22:19:09.0 


+56:04:45.0 


Pal 


148 


22305+5803 


22:32:24.3 


+58:18:58.2 


Pal 


151 


22506+5944 


22:52:38.6 


+60:00:56.0 


Pal 


160 


23385+6053 


23:40:53.3 


+61:10:19.1 


Pal 



" Source running number from Molinari et al. ll996r 
^ Imaged only in Ks . 



The first approach we tried to follow is the PSF-fitting pho- 
tometry which should be less prone to these problems. We chose 
a sub-sample of test fields with diff'erent levels of stellar crowd- 
ing. An important aspect in this procedure is the modelling of 
the PSF; we made several trials selecting a variable number of 
point-like sources (from 3 to 30) of diff'erent brightness levels 
and diff'erent position in the field. We find that the resulting 
PSF model is not particularly sensitive to the choice of numbers 
and/or brightness of the stars; on the other hand the results are 
quite diff'erent depending on the mean stellar density of the field. 
The photometry was carried out using the ALLSTAR task, par- 
ticularly suited for crowded fields. However, we also tested the 
other two tasks (PEAK and NSTAR) and they produce compara- 
ble results for most of the sources. We note, however, that espe- 
cially in the most crowded areas the subtraction of the PSF-fitted 
sources from the image introduces two spurious eff'ects: an unac- 
ceptably high level of residuals with brightness levels well above 
the detection threshold used and a great number of negative holes 
which show that the psf-fit includes a portion of background in 
the extimation of the source flux, therefore overestimating its 
value. This is due both to the limited accuracy of the PSF model 
that can be built in very crowded fields, where faint neighboring 
stars can enter in the area where the PSF model is estimated, and 
to the presence of a significant and variable background which 
is quite common and expected in the Galactic Plane. A similar 
conclusion was reached by Hillenbrand & Carpenter (120001) in 
their study of the inner Orion Nebula Cluster. 

The second approach we followed is standard aperture pho- 
tometry. The choice of the radii for the aperture and for the back- 
ground annuli is of course extremely important. The optimum 
aperture should not be too large to include nearby sources and 
not too small to significantly cut the PSF and severely under- 
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estimate the flux. We did several attempts on one of the most 
crowded fields (Mol30, observed at NTT) with three diff'erent 
aperture radii equal to the PSF FWHM (typically O'/l at NTT 
and I'/ 4 at Palomar in K^), twice and thrice this value. For each 
photometry run we analyze the sources' flux distribution and, 
as expected, the median flux is found to increase with increas- 
ing aperture radius. Increasing the aperture from one to two PSF 
FWHMs raises the median source flux by an amount compatible 
with the inclusion of the first ring of an Airy diffraction pattern; 
instead, an aperture radius equal to three times the PSF FWHM 
produces a ffux increase much too large compared to the addi- 
tional fraction of the Airy profile entering the aperture, therefore 
being caused by the inclusion of nearby sources. We adopted an 
aperture radius equal to the PSF FWHM to minimize neighbour 
contamination, and then applied an aperture correction factor for 
the fraction of the PSF cut out of the aperture; this was estimated 
via multi-aperture photometry (starting from a 1 FWHM size) on 
relatively isolated stars in the target fields. 

A further eff'ect to be corrected for, also given the crowd- 
ing of our fields, is the possible contamination arising from the 
tails of the brightness profiles of neighbouring stars. To quan- 
tify this contamination we created a grid of simulations with two 
symmetric Gaussians with a wide variety of peak contrast and 
at diff'erent reciprocal distances. We computed the fraction of 
the Gaussian profile of the neighbouring source that falls into 
the photometry aperture centered on the main source, and hence 
generated a matrix of photometry corrections for diff'erent source 
distances and peak contrasts. At this point we run throught the 
magnitude file produced by the aperture photometry task and for 
each source we apply a magnitude correction depending on the 
presence, distance and contrast ratios with other neighbouring 
stars. 

In spite of the various issues discussed above, the photome- 
try obtained with the two methods are in a good agreement with 
each other, except at faint magnitudes. For these faint objects 
we always find that the PSF photometry tends to produce lower 
magnitudes (hence stronger sources) than the aperture photome- 
try; this eff'ect is easily understood given our finding (see above) 
that the subtraction of PSF-fitted sources always leaves negative 
holes in the residual image, and this eff'ect is much more im- 
portant for faint stars. We thus decided to adopt the magnitudes 
determined from aperture photometry. 

For each target field we estimated the limiting magnitude 
(LM) using artificial star experiments. The fields were populated 
using the IRAF task ADDSTAR with 400 fake stars with mag- 
nitudes distributed in bins of 0.25 mag between values of 15 and 
21; the percentage of recovered stars as a function of magnitude 
gives an estimate of the completeness level of our photometry. 
The star recovery percentage does not decrease monotonically 
with increasing magnitude because fake stars can be placed also 
very close to bright real stars and then go undetected by the find- 
ing algorithm. However, we find that the limit of 85-90% recov- 
ery fraction is reached on average around J=18.7, H=17.7 and 
K,=17.4 for NTT images, and J=18.0, H=17.3 and K,=16.6 for 
Palomar images. We find that the typical photometric uncertainty 
is below O.lmag close to the limiting magnitude. 

To verify the integrity of our photometry we compared our 
magnitudes with those extracted from 2MASS point source cat- 
alog for all the fields of our sample. Considering the diff'erent 
spatial resolutions between 2MASS and the telescopes used for 
our observations, this comparison was limited to those 2MASS 
point-like sources associated with a single source in the Palomar 
or NTT images. The median diff'erences with respect to 2MASS 
for the various fields are of the order of -0.1, -0.2 and -0.3 mag 



for J, H, and bands, respectively. Within each field, the scatter 
around these median values is ~ 0.1 mag in all three bands, con- 
firming the internal consistency of our photometry. Noticeable 
departures (-0.5 mag) of the median diff'erence with 2MASS 
from the above values are observed for the field of source Moll 1 
(Palomar), and for sources Moll03, Moll09 and Moll 10 (NTT). 
However, the latter sources were observed on the same night, 
which our log registered as not good due to sky variations which 
are not tracked by night-averaged zero-points. We stress again, 
however, that these are systematic diff'erences with respect to 
2MASS in this limited number of cases; the r.m.s. scatter around 
these median diff'erences are ~ 0.1 mag in all bands and this 
should give confidence that the internal consistency of the pho- 
tometry in each field is preserved. We then decided to rescale 
our photometry to the 2MASS photometric system to remove 
these systematic eff'ects. The (J-H) and (H-K) color diff'erences 
between 2MASS and our photometry are not correlated with the 
magnitude, so that no magnitude-dependent color eff'ect is intro- 
duced in this rescaling. 

3. Results 

3.1. Cluster Identification 

The identification of a cluster results from the analysis of stellar 
density in the field. Since our target fields are sites of massive 
star formation associated with local peaks of dust column den- 
sities and hence of visual extinction, the Ks images are clearly 
more suited for this type of analysis. 

Stellar density maps were built for each field counting stars 
in a running boxcar of size equal to 20''. The box size was deter- 
mined empirically to enhance the statistical significance of local 
stellar density peaks and to maximize the ability to detect the 
clusters. Larger boxes tend to smear the cluster into the back- 
ground stellar density field decreasing the statistical significance 
of the peak, which may lead to non-detection of a clearly evi- 
dent cluster, particularly in the rich inner Galaxy fields (this hap- 
pens, e.g., for source Moll03, see Fig. |A.18] in the Appendix). 
Smaller boxes produce noisy density maps where the number of 
sources in each bin starts to be comparable to the ffuctuations of 
the background density field due either to intrinsic variations of 
the field star density or to variable extinction from diff'use fore- 
ground ISM in the Galactic Plane (where all of our sources are 
located). For most of our objects in the outer Galaxy this analy- 
sis is just used to locate the position of the peak stellar density, 
since the clusters are obvious already from the visual inspec- 
tion (Mol3 to Mol28, and Moll43 to Moll51, see Appendix). 
For the rest of the fields the density maps are used to ascer- 
tain the presence of a cluster; especially toward the inner Galaxy 
the density maps tend to show more than one peak at compa- 
rable levels. It is important to remember, however, that this is 
a search for stellar clusters toward regions where indications of 
active star formation are already available, and this information 
can be used. In particular, the coincidence of these peaks with 
cold dust clumps as traced by intense submillimeter and mil- 
limeter emission (Beltran et al. 2006 , Molinari et al. "2008a) is 
critical to consider the density peak as a real feature associated 
with the star formation region. Casual association is ruled out by 
the high number of positive associations (see Table |2]). 

As a further confirmation step for the positive detection of 
a cluster we build radial stellar density profiles where stars 
are counted in annuli of increasing internal radius and constant 
width and then ratioed to the area of the annuli (Testi, Palla & 
Natta J998j); uncertainties are assigned assuming Poisson statis- 
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Fig. 1. Stellar density (in stars/pc^), for Mol28, as a function of 
the radial distance (in parsecs) from cluster center. Error-bars are 
computed as the Poissonian fluctuations of source counts in each 
bin. 



tics on the number of stars in each annulus. We then assign a 
positive cluster identification if the radial profile exhibits at least 
two annuli with values above the background. In order to re- 
fine the location of the density peak, we repeat the radial density 
profile analysis starting from several locations within W of the 
peak derived from the density maps; the location which maxi- 
mizes the overall statistical significance of the annuli is then as- 
signed as the cluster center. Figure [T] shows the typical footprint 
of a cluster, where the stellar density is plotted as a function of 
distance r from the start location; the density has a maximum at 
r=0 and decreases until it reaches a constant value which is the 
average background/foreground stellar density. 

There are two exceptions in this analysis. The first is for 
source Moll 60. The K^-band image shows clear stellar density 
enhancement in a semi-circular annulus which surrounds the 
northern side of the dense millimeter core (see Figure |A.26| in 
the Appendix), which appears devoid of stars. This stellar den- 
sity enhancement is coincident with the emission patterns visi- 
ble in the mid-IR (Molinari et al. 2008b ) so it is clearly a stel- 
lar population associated with the star forming region. Since the 
millimeter peak is at the center of symmetry of the semi-circular 
stellar distribution, we will consider this as the center of the clus- 
ter. This is just for completeness of reporting, since we cannot 
say if the low density of stars at the millimeter peak is an ef- 
fect of extreme visual extinction or reflects an intrinsic paucity 
of NIR- visible forming stars, as the proposed extreme youth of 
the massive YSO accreting in its depth would seem to suggest 
(Molinari et al. 2008b). 

The second exception is for source Mol8. The stellar density 
analysis shows two peaks which are coincident with two distinct 
dust cores (see Fig. A.2); we then assume the presence of two 
distinct clusters, rather than a sub-clustering feature within the 
same cluster. The radial density profile analysis cannot be used 
here, so we fit elliptical gaussians to the peaks in the density 
maps, allowing for an underlying constant level representing the 
background stellar density. The resulting cluster richness is ob- 
tained integrating the fitted gaussian, and the cluster radius is 
taken equal to the fitted FWHM (the fitted gaussians are nearly 
circular). 



Always following Testi et al. (119981) , we determine the rich- 
ness indicator of the cluster 1^ by integrating the background- 
subtracted density profile; the cluster radius is taken as the radial 
distance from start location where the density profile reaches 
a constant value. This richness indicator is a very convenient 
figure to use in these cases where no detailed information is 
available for each single star in the region and the member- 
ship of the cluster cannot then be established for each single 
star. These values are reported in Col. 3 of Table |2] for all 
fields where a cluster has been clearly revealed. Col. 1 gives 
the target name (cf. Table [T]); it's kinematic distance is listed 
in Col. 2. The parameter Nobs (Col. 4) is the number of clus- 
ter members derived (see Sect. |4.1| below) from the integration 
of the background- subtracted Kg Luminosity Function (here- 
after KLF, see Sect. |4.1| ). Also reported in Col. 8 is the mass 
of the hosting molecular clump; this is derived from the cold 
dust emission as reported in Molinari et al. 2008a, 2000), in- 
tegrated over the entire spatial extent of the cluster; conversion 
into masses is done under opticall y thin a ssumption assuming 
T=30 K, 13 = 1.5 (Molinari et al. 2008a) and a mass opacity 
f<230GHz = O.OOScm^g'^ which incorporates a gas/dust weight 
ratio of 100 (Preibisch et al. 119931 . The IRAS source bolo- 
metric luminosity. Col. 9, is taken from Molinari et al. (1 9961 
I20QQ1 120021 l20Qgal) : in Col. 10 we fist the Ay on the peak clus- 
ter position as estimated from submm observations (Molinari 
et al. 2008a, 2000). In Col. 11 and Col. 12 the coordinates of 
centers of the identified clusters are reported. Columns 6 and 7 
contain parameters that will be described later in the text (see 
Sect.[T2]). 

Following the described procedure, a cluster is detected 
within V of the IRAS position for 22 out of the 26 observed 
fields (85% detection rate). In two cases (Mol38 and Mol59) the 
stellar density map does not show a clear peak above the fluctu- 
ations of the field stellar density. For Mol 98 the radial density 
profile only shows one annulus above the background, so they 
fail the criterion that the stellar density enhancement should be 
significantly resolved above the background in two annuli. In 
one case (Mol30) several stellar density peaks are found in prox- 
imity of the IRAS source, but the lack of information of submil- 
limeter/millimeter continuum prevents any firm conclusion. 

Figure [2] shows the run of 1^ as a function of peak Ay and 
suggests that with larger dust extinction it may be more difficult, 
or less likely, to detect a cluster at 2.2 jim. 

Our detection rate is quite high and tells us that young stellar 
clusters in sites of intermediate and massive star formation are 
essentially ubiquitous. While this evidence was established for 
relatively old Pre-MS systems like Herbig Ae/Be stars (Testi et 
al. 1999), we hereby verify that this is the case also in much 
younger systems where the most massive stars may even be in a 
pre-Hot Core stage (Molinari et al. l2008aj) . 

Our detection rate is higher compared to other similar 
searches of stellar clusters toward high-mass YSOs. For example 
Kumar et al. ( 2006 ) use the 2MASS archive and report a rate of 
25% (rising to 60% neglecting the inner Galaxy regions) toward 
a larger sample which also includes the sources of this work; in 
particular we detect all clusters also detected by Kumar et al. and 
in addition we reveal clusters toward 13 objects for which Kumar 
et al. report no detection. The reason for this discrepancy may 
be due to the fact that we obtained dedicated observations while 
Kumar et al. used data from the 2MASS archive; the diffraction- 
limited spatial resolution of our data is between a factor of 4 and 
a factor of 10 better with respect to 2MASS, and this certainly 
facilitates cluster detection especially in particularly crowded ar- 
eas like the inner Galactic plane. To test this hypothesis we de- 
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Sou. 






^obs 




Pre-MS 




^bol 


Ay peak 


Cluster Center 


Mol 


kpc 






pc 


% 
CC 


% 

CM 


Mo 


lO^Lo 


mag 


Qf(J2000) 


(5(J2000) 


3 


2.17 


78 


78 


1.7 


34 


99 


910 


12.4 


18 


00:44:57.4 


+55:47:20.0 


8A 


11.5 


25 


30 


1.3 


37 


27 


1650 


57.0 


18 


05:17:13.8 


+39:22:29.7 


8B 


11.5 


27 


24 


1.3 


9 


4 


1780 


5.5 


8 


05:17:12.0 


+39:21:51.8 


9 


6.2 


7 


7 


0.6 


16 


_<? 


~ 


24 


~ 


05:20:16.9 


+36:37:22.0 


11 


2.1 


51 


48 


0.5 


12 


41 


360 


4.6 


35 


05:37:47.7 


+31:59:24.0 


12 


1.6 


12 


13 


0.3 





30 


72 


1.6 


46 


05:40:24.4 


+23:51:54.8 


15 


1.5 


64 


61 
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^ Kinetic distance using the rotation curve from Brand & Blitz' 1993' 
^ Stellar density analysis inconclusive due to extreme crowdedness of this field. 
^ Stellar density reveales no peaks close to the IRAS position or the submm peak. 
^ Only observed in Kg. 

^ Detection refused due to extreme field complexity (see text). 

^ Detection refused because only 1 annulus in the radial density profile is above background (see text). 

^ No extinction estimate is available due to lack of submm information to evaluate de-reddening correction. 

^ Extinction estimate is available from single-pointing submillimeter data (Molinari et al. 2000J but not from maps, so that a 

reliable clump mass estimate is not possible. 




Fig. 2. Cluster richness indicator 1^ as a function of Ay on cluster 
center; for a few detected cluster we do not have an estimate of 

Ay). 



graded the NTT image of Mol 103, also considered in Kumar 
et al., to the 2MASS resolution; extraction and photometry were 
performed as outlined above but the search for a cluster based on 
the stellar radial density profiles revealed no cluster. Besides, the 
estimated number of members (corrected for the contribution of 
fore/background stars) for 7 out of the 10 clusters detected both 
by us and by Kumar et al. is at least a factor of two less in the 
latter study. 



Kumar & Grave (12008b conducted a similar study on a large 
sample of high-mass YSOs, that include a certain number of 
our sources, using this time data from the GLIMPSE survey 
(Benjamin et al. 2003 ). In this work they detect no significant 
cluster around any targets in a sample of 509 objects. As the au- 
thors say in the paper, however, GLIMPSE data are sensitive to 
2-4 M© pre-main sequence stars at the distance of 3 kpc. Based 
on color-magnitude analysis (see later below) our mass sensitiv- 
ity is of the order of 1 M© at a distance of 3.6 kpc and -0.6 M© 
at a distance of 2.1 Kpc. Probing longer wavelengths, GLIMPSE 
is likely to be more sensitive to younger sources compared to 
the classical J, H, K range which also samples relatively older 
pre-MS objects. The combination of sampling higher-mass (and 
hence more rare stars due to the shape of the IMF) and relatively 
younger stars (which, as indeed our analysis finds, may not be 
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the majority in a young cluster) may plausibly be the reason of 
the negative cluster detection results in Kumar & Grave. 
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Fig. 3. Distribution of the cluster's radii in parsecs (full line), and 
of the cluster richness indicator 1^ in number of stars (dashed 
lines); the median values for the two distributions are 0.6 pc and 
37 stars, respectively 

The distribution of the radii of the detected clusters is shown 
with the full line in Fig.[3j the median value is 0.7 pc. The dashed 
histogram (which refers to the upper X-axis) shows the distribu- 
tion of the cluster richness indicator 1^, with a median number of 
stars of 27. We note that the value of 1^ for many of our clusters 
is less than the limit of 35 suggested by Lada & Lada ( 2003 ) to 
be a bona fide cluster. This definition stems from the argument 
that a less rich agglomerate may not survive the formation pro- 
cess as an entity. Our interest, however, is to investigate the spa- 
tial properties of the young stellar population in a star forming 
region at the time of active formation, without worrying about 
its possible persistence as a cluster at the end of the formation 
phase. However, we prefer not to introduce a new term to iden- 
tify the structures we see and still use the term cluster, although 
in a ''weaker'' sense compared to Lada & Lada. 

3.2. Properties of Identified Clusters 

We will first derive qualitative indications concerning the nature 
of the identified clusters using simple diagnostic tools like color- 
color and color-magnitude diagrams. These diagrams have been 
drawn for all detected clusters and are available in electronic 
form; we here illustrate the particular case for Mol28. 



cal of young pre-MS objects with an intrinsic IR excess aris- 
ing from warm circumstellar dust distributed in disks (Lada & 
Adams 1992). The set of dotted curves represents the locus of 
two-component black-bodies with temperatures as indicated at 
the start and end of each dotted line; along each curve the relative 
contribution of the the two blackbodies is varied. These curves 
mimic the eff'ect of a temperature stratification in the dusty cir- 
cumstellar envelopes, and the presence of sources in the area 
covered by these curves is an indication of the presence of warm 
circumstellar dust. 
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Fig. 4. [J-H]v^[H-Ks] diagram for Mol28. [J-H] is obviously a 
lower limit for sources not detected in J. The continuous curve at 
the bottom-left represents the Main Sequence, while the dashed 
grey lines represent the eff'ect of reddening (Rieke & Lebofsky 
1985 ) for variable amounts of extinction as indicted along the 
lines. The dashed-dotted black line is the Black-Body curve, and 
the dotted curves are two-component Black-Body with varying 
relative contribution (respectively, from the inner to the outer 
curve, 3000-1500K, 3000-lOOOK, 3000-900K and 3000-500K). 

A straightforward indication of the youth of the cluster may 
be off'ered by the fraction of sources which are not compatible 
with reddened MS stars, i.e. those with IR excess. The number 
of stars with IR excess is normalized to the total number of stars 
revealed in the cluster area corrected for the expected number of 
fore/background stars estimated from the areas surrounding the 
cluster (but still in the same imaged field). To be conservative 
we extend the region of the MS by 0.2 magnitudes to the right 
corresponding to about 2cr uncertainty on measured magnitudes. 
This ratio is reported as a percentage value in Col. 6 of Table [2]). 



3.2.1. Color-Color Analysis 

Fig. |4] shows the [J-H]v^[H-Ks] diagram for all sources detected 
within a distance equal to R^/ centered on the stellar density 
peak. The full circles represent all sources detected in all three 
bands, the arrows represent sources with lower limits (in magni- 
tude) in the J band. The plot shows more stars than the 1^ value 
reported in Table [2]because we also include the fore/background 
stars that cannot be individually identified against the true clus- 
ter members. A significant fraction of the sources have colors 
compatible with main- sequence stars with a variable amount of 
extinction reddening (computed adopting the Rieke & Lebofsky 
((1985i) extinction curve), but many sources show colors typi- 



3.2.2. Color-Magnitude Analysis 

Additional evolutionary indications for the detected clusters may 
be derived from the Ks-[H-Ks] diagram, reported for Mol28 in 
Fig.|5] Compared to the main sequence (the leftmost almost ver- 
tical curve in the figure) a significant fraction of the sources are 
on its right, where the evolutionary tracks for Pre-MS sources 
(Palla & Stabler 119991 ) can also be found, and could be there- 
fore interpreted as very young pre-MS objects. The distribution 
of sources in the diagram spans a much larger region that the one 
covered by the Pre-MS isochrones, and this is due to the com- 
bined eff'ect of extinction reddening and IR excess. Extinction 
eff'ects can be appreciated looking at the dotted lines original- 
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ing on the main sequence and extending toward bottom-right 
for increasing values of Ay; on the other hand the presence of 
a warm dusty circumstellar envelope implies an increase in ab- 
solute emission and in SED steepness, therefore shifting a pure 
photosphere toward top-right of the diagram (as shown by the 
arrow labeled 'IREX' in Fig. 5]). Similar to the color-color anal- 
ysis, it is impossible to try and estimate the age of individual stel- 
lar sources based on their location on the pre-MS isochrones, be- 
cause we do not know the amount of Ay by which we should de- 
redden each object. We follow a conservative approach by dered- 
dening each object using half of the exctinction estimated for 
each location from millimeter maps; this corresponds to putting 
each object midway through the clump. 
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Fig. 5. KsV^[H-Ks] diagram for Mol28. The leftmost curve repre- 
sents the Main Sequence, while the dashed lines represent the ef- 
fect of reddening for variable amounts of extinction. Isochrones 
from Palla & Stabler (1999 ) are also indicated with full lines for 
different Pre-MS ages. The arrow labeled IREX indicates the di- 
rection of change due to IR-Excess (see Sect. |4.2| ). [H-K] is obvi- 
ously a lower limit for those sources not detected in H. Symbols 
in grey color indicate sources with IR excess as determined from 
the color-color diagrams (see fig.|4]). 

A further correction is to remove the IR excess for those 
sources in which this is apparent in the color-color diagram 
(see fig. |4]), as estimated using the formulation suggested by 
Hillenbrand & Carpenter (I2QQ0I) , and that will also be used 
later in this work (see Sect. |4.2| ). The fraction of pre-MS stars 
over the total in each cluster area will still be contaminated by 
fore/background stars; to estimate this contamination we choose 
an off'-cluster area in the same imaged field and we simply com- 
pute the ratio of sources with pre-MS colors over the total (in 
these off'-cluster regions there is no significant reddening to cor- 
rect for). Col. 7 of Table [2]reports for each cluster the fraction of 
stars (detected in the cluster area in all three bands) situated more 
than 0.2 mag to the right of the MS after the various corrections 
have been applied. 

4. Initial Mass Functions and Star Formation 
Histories 

As it is apparent from the qualitative analysis presented in the 
previous paragraphs, the diagnostic power of our observations is 
limited because we do not know which objects in the cluster area 



are real cluster members and we do not know the exact amount 
of dust extinction (originating within the hosting clump) and 
IR excess (originating in the immediate circumstellar environ- 
ment) pertaining to each source. Lacking the detailed knowledge 
on individual stars in the clusters, fundamental quantities like 
the Initial Mass Function (IMF) and the Star Formation History 
(SFH) cannot be directly derived from, e.g., the Kg luminosity 
function (KLF). We are then forced to obtain these using statis- 
tical simulations of clusters based on diff'erent input parameters 
and performing a statistical comparison between synthetic and 
observed KLFs and HKCFs. 

We will first derive the observed KLFs from the observa- 
tions; we will then illustrate in detail the model used for the clus- 
ter simulations, exploring the sensitivity of the results to a wide 
range of input parameters; finally, modeled and observed KLFs 
will be compared to infer statistical indications for the IMF and 
SFH for our clusters. 



4.1. Observed Luminosity Functions 

The KLF for each cluster is obtained simply counting all de- 
tected sources within th e clu ster area as identified from the clus- 
ter density profile (see |3.1| ). Similar to the other diagnistic tools 
( |3.2.1| and |3.2.2| ), the KLF will be contaminated by field stars 
that cannot be individually identified. To account for the field 
star contamination in a statistical way we subtract from the KLF 
built on the cluster area, the KLF built in a region outside the 
cluster area but still in the same imaged field, after normalis- 
ing for the diff'erent areas. Regions where the field star KLF is 
built have a lower extinction respect to the cluster one, so the 
background contribution to the cluster KLF is likely to be over- 
estimated. Field- subtracted KLFs for all clusters are presented 
in Appendix [5] and are also available onlin^ 

The integral of the KLF gives an independent estimate of the 
number of cluster members; these values are reported as Nobs 
in Table |2j their agreement with richness indicator 1^ confirms 
the consistency of the analysis. All KLFs show a dominant peak 
which always lies close to the completeness limit, showing that 
our observations are not sensitive enough to trace the low-mass 
stellar component of our clusters. Many of the KLFs present a 
separate small peak at low magnitudes (one or two sources at 
most, on average). Can this be due to confusion arising from 
source crowding and insufficient spatial resolution ? We studied 
for each cluster the distribution of distances of each star from its 
nearest neighbour and we find that there are essentially two types 
of distributions, reported in Fig. [6] In the first type (full line in 
figure) the distribution has a peak corresponding to an inter- star 
distance significantly higher than the value corresponding to half 
the PSF FWHM (the full vertical line); in this case the suggestion 
is that all cluster members have been resolved from their neigh- 
bour. In the second type (dashed line in the figure) the distribu- 
tion has its peak very close to half the PSF's FWHM (the dashed 
vertical line), indicating that source blending should be certainly 
considered possible. We verified that all clusters exhibiting a dis- 
tance distribution of the second type do show a second faint peak 
at high brightness in their KLFs, therefore confirming that this 
feature is an artifact of the relatively low spatial resolution which 
in some cases is not sufficient to resolve all cluster members. 



at http ://galatea.if si-roma.inaf.it/faustini/KLF/ 
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repeated until the number of synthetic detectable stars equal the 
value of Ic as determined in our observations; at this point the 
cluster generation process is complete. 

Since the simulation is based on Monte-Carlo extraction for 
stellar mass, age and position in the cluster, each independent 
run for a fixed set of input parameters can in principle result in 
very diff'erent outputs in terms of cluster luminosity, total stellar 
mass, maximum stellar mass and synthetic KLF. To provide sta- 
tistical significance, the model is run 200 times for any given set 
of input parameters, and the median KLF is later adopted for the 
comparison with the observed one. Clearly, the predictive power 
of this simulation model resides in its capability to characterize 
the cluster properties for any given parameter set. In other words, 
the distribution of the resulting quantities should not be uniform 
but peaked around characteristic values. We will come back to 
this in section Sect. 14.2.3] 



Fig. 6. Distribution of identified sources as function of nearest- 
neighbour distance (D) for two of examined fields (Mol28 
dashed line and Mol 103 full line). 

4.2. Synthetic KLF. Synthetic Cluster Generator: a Near-IR 
Cluster Simulator 

As already mentioned, from our data alone we cannot derive 
masses and ages. We thus developed a model to create statis- 
tically significant cluster simulations obtained for diff'erent as- 
sumptions of IMF and SFH (source ages and their distribu- 
tion), and compare the synthetic KLFs with the observed field- 
subtracted KLFs. This model we call Synthetic cluster Generator 
(SCG). 

4.2.1. SCG: Model Description 

A cluster is created by adding stars whose masses and ages are 
assigned via a Monte-Carlo extraction according to the chosen 
IMF an d SFH; the pre-MS evolutionary tracks of Palla & Stabler 
(119991) are then used to convert them into J, H and Kg mag- 
nitudes. The 3D distribution of stars is obtained by randomly 
choosing for each star a set of xj,z coordinates using the ob- 
served stellar density profile (see Sect. |3.1| ), approximated as a 
radially symmetric Gaussian, as weight- function; this is needed 
to assign, using submm continuum images, the proper column 
of cold dust to extinguish the near-IR radiation. Other analytical 
functions could have been used, e.g. a King profile, but the statis- 
tics of our clusters are not sufficiently high to try and explore the 
eff'ect of diff'erent radial profile assumptions. 

To convert the submm flux into dust column density we used 
the dust temperature and emissivity exponent JS as determined in 
Molinari et al. (2000); mean values from the latter work were 
adopted for those fields not covered as part of that work. 

To properly simulate the pre-MS stars, we also need to in- 
clude the effect of an IR excess due to warm dust in the circum- 
stellar envelopes and disks. We used the distribution (modeled as 
a Gaussian) of [H-Kslg;^ color excesses as measured for a sample 
of Pre-MS stars in Taurus, as used by Hillenbrand & Carpenter 
(120001) , as a weight-function to randomly assign a [H-KsJ^x to 
each simulated star in our model; the KsV^ [H-Kslg;^ relationship 
adopted in the above mentioned work was then used to derive 
the H and Kg excess-corrected magnitudes. The magnitude 
of the synthetic star was then compared with the limiting magni- 
tude typical of the cluster being simulated to determine if the star 
could have been detected in our observations. This procedure is 



4.2.2. SCG: Input Assumptions 

We tested three different assumptions for star formation histories 
in our cluster simulations. The first choice is to assume that stars 
in the cluster formed in a single burst- like event (hereafter SB) 
some ti years ago. The explored range in the simulations is 10^ < 
ti < 10^ yrs. The second choice is to have the formation of stars 
proceed at a constant rate (hereafter CR) from a time ti years ago 
to a time t2 years ago. The ranges explored in the simulations are 
10"^ < ti < 10^ yrs and 10^ <t2< 10^ yrs, where always ti > t2. 
The third possibility we explored is a variation of the previous 
one, where the star formation rate is not constant but varies with 
time in a Gaussian fashion (hereafter GR). Within the boundaries 
for start and end of the star formation process, ti and t2 varied 
as above, we also varied both the time t^ for the peak of the 
Gaussian in the range 10^- < t^ < 10^-^; the Logio(cr) of the 
Gaussian-like SFH was allowed to assume the two values 0.1 
and 0.5. 

As for the IMFs, we allowed three different choices from 
Kroupa et al. ( 1993 ), Scalo ( 1998 ) and Salpeter ( 1955 ), with the 
latter modified introducing a different slope for M<1 M© that 
coincides with that of the Scalo (1998 ) IMF; the three IMFs 
will be coded as IMFl, IMF2, and IMF3 respectively. The IMF 
from Kroupa et al. favors the low-mass end of the distribution, 
while the classical Salpeter IMF is flatter at low mass but heav- 
ier at intermediate and high masses (above 1 M©). The IMF 
from Scalo is kind of intermediate between the two, resembling 
Salpeter 's one below 1 M© and above 10 M©, and Kroupa' s for 
1 M©<M<10M©. 

4.2.3. SCG: Predictive Power 

In order to verify our model's predictive power, we ran a set 
of 200 simulations for a cluster with a Salpeter IMF and a con- 
stant star formation rate with ti = 10^ yrs and 12=10"^ yrs. Figure]?] 
shows the distribution of the predicted number of stars and the 
total luminosity for the 200 simulations. The number of cluster 
members shows very little variations, as expected since the num- 
ber of detectable stars is the parameter we use to stop the simu- 
lation; on the other hand the distribution of the total luminosity 
is not particularly peaked, as the central 3 bins containing about 
60% of the simulations span almost two decades in luminosity. 

On the other hand the distributions for the total cluster stel- 
lar mass, and for the mass of the most massive member (see 
Fig. [8]) are rather peaked and highlight a relatively higher pre- 
dictive power of the model for these two quantities. It is to be 
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Fig. 7. Distribution of the predicted number of cluster members 
(full line) and total luminosity (dashed line) for 200 SCG runs 
for Moll 60 with a Salpeter IMF and a constant star formation 
rate with ti = 10^ yrs and 12=10"^ yrs. 



noted that the distributions are rather skewed, suggesting that 
neither the mean nor the median are particularly suited to char- 
acterize the peak of the distribution. Indeed we have computed 
the values that these quantities assume at the peak of their distri- 
butions, to have a more representative value for the masses and 
use them in the following arguments. 
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Fig. 8. Distribution of the predicted total stellar mass (full line) 
and mass for the most massive star (dashed line) in a cluster for 
200 SCG runs for M0II6O (same inputs as in Fig.^. 

Concerning the reproducibility of the KLF, for each of the 
200 runs the resulting KLF was fitted with a Gaussian func- 
tion and the center, peak and cr were determined. Fig. [9] reports 
the distribution of these three parameters for the 200 runs and 
shows that the all of them are remarkably peaked and symmet- 
ric. The formal r.m.s. spread for the three quantities, estimated 
via a Gaussian fit to the distributions in the figure, is 0.3 mag 
for the KLF center, ^ 12% for the KLF peak (about 1.2 sources 
out of a mean KLF peak of 10), and ^0.25 mag for the KLF 
FWHM. 



Fig. 9. Distribution of the predicted center magnitude (full line - 
bottom X-axis scale), width (dotted line - bottom X-axis scale) 
and peak value (dashed line - top X-axis scale) of the predicted 
Gaussian-fitted KLFs for 200 SCG runs for Moll 60 (same inputs 
as in Fig. [7]). 



We have made a similar analysis for HKCF (H-Kg color 



10 shows the distribution of 



function; see Sect. 3.2.2). Figure 
Gaussian function centers, peaks and cr's for HKCFs obtained 
for the same 200 runs used previously for the KLFs. Gaussian 
fits to the three distributions in the figure give an r.m.s. that is 
^0.15 mag for the HKCF center and ^ 0.14 mag for the HKCF 
FWHM, while the 'peak' distribution is flatter and has an r.m.s. 
value of 21% for the HKCF peak (about 3.2 sources out of a 
mean HKCF peak of 15). It is worthwhile to stress that since the 
position which is assigned to each simulated star in the cluster 
is difl'erent in each of the 200 runs of the model (for any given 
set of input parameters), the scatter in the properties of the syn- 
thetic KLFs and HKCFs also statistically tends to account for 
the efl'ects of extinction variations in the cluster's hosting clump, 
which may in principle be relevant in such heavily embedded 
systems (see Tab.|2]). 
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Fig. 10. Distribution of the predicted center color (full line - 
bottom X-axis scale), width (dotted line - bottom X-axis scale) 
and peak value (dashed line - top X-axis scale) of the predicted 
Gaussian-fitted HKCFs for 200 SCG runs for Moll 60. 
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We conclude that the model results, for a given set of input 
parameters, have a good reproducibility, except concerning the 
total luminosity. The model therefore has a quite fine predictive 
power concerning the median properties of a synthetic cluster. 
Indeed, the spread of KLF center magnitudes is less than the bin 
amplitude used in building the KLFs for the simulations (and 
that will be used for the rest of the work); the median synthetic 
KLF therefore provides a good representation of the cluster lu- 
minosity distribution. 

In conclusion, 200 simulation runs for each combination of 
input parameters (IMF and SFH) seem enough to achieve a suf- 
ficient statistical significance for the properties of the synthetic 
observables (KLFs and HKCFs). Although the distributions for 
the KLFs' (HKCFs') parameters seem rather symmetrical, we 
will adopt the median KLF (HKCF) of the 200 runs as a more re- 
liable characterization for that particular parameters' set. Using 
the mean KLF (HKCF) for the comparison does not significantly 
alter the results. 

4.2.4. Exploring the SCG parameter space: Cluster 
Parameters 

After verifying the robustness of model results for independent 
runs for the same input parameters, we now want to verify the 
sensitivity of the model results against changes of these param- 
eters. We will first concentrate on simulated cluster physical pa- 
rameters (number of cluster members, total luminosity, stellar 
mass distribution), and in the next paragraph we will examine 
how the KLFs and the appearance of the color-magnitude dia- 
grams, which are the main observables we will use in our analy- 
sis, behave in this respect. 



ber of star extractions and hence a lower relative total mass at 
the end of the simulation. The age-trend of M^^Max is instead the 
same (only shifted toward higher masses) because the probabil- 
ity of extracting a massive star is the same for all ages and is 
only function of the chosen IMF. 

Total Stellar Luminosities and Massive Object Luminosities - 

The total stellar luminosity, like the luminosity of the most mas- 
sive star (hi, Max) r exhibits the same behaviour as M^max- This is 
easily understood given the steep power-law dependence of the 
stellar luminosity on mass, and confirms that the total luminos- 
ity (hrot) will be largely dominated by the most massive stellar 
object in the cluster: Lxot ^ ^i^Max- Of great interest is the ra- 
tio between h^Max and hrot', for the large majority of clusters its 
value varies between 0.6 and 0.8. This is further confirmation 
that global properties of our clusters are dominated by the most 
massive source. This ratio does not present any particular depen- 
dence on the value of M^max, or of N^^^^ ; only for the most pop- 
ulated clusters (clusters with of the order of a hundred members, 
such as Moll 03, where the contribution of a great number of 
low-mass sources becomes important, do we find a lower value 
for this ratio. 

4.2.5. Exploring the SCG parameter space: KLF variations 

We will now briefly analyze the diagnostic power of the KLF and 
the HKCF against changes in IMF and SFH choices. Figure [TT 
shows the KLFs predicted for source Mol3 adopting the same 
SFH parameters (as indicated in the figure) and using the three 
diflferent IMF choices. 



Number of stars l^stars - As a general rule the older the cluster 
is allowed to be, irrespective of the detailed SFH adopted, the 
higher is the number of produced stars. This is easily understood 
since the SCG cluster building is stopped when the number of 
the Ks -detectable stars equals the number of observed objects; 
if a cluster is old the stars will be intrinsically fainter due to the 
shape of Pre-MS tracks, and it will statistically be less likely 
to extract stars bright enough to be detectable. N^^^^ does not 
significantly depend on the IMF choice as long as ti < 10^ yrs, 
while for older systems IMFl (Kroupa et al. 1993 ) will produce 
nearly twice as many stars as IMF3 (Salpeter J955J with IMF2 
(Scalo [T998l in between. 




Stellar Masses - Likewise, the total stellar mass and the mass 
of the most massive star will be higher the older the cluster is 
allowed to be. If an IMFl cluster is a very old SB or a CR with 
ti=10^ yrs and t2=10 yrs for example, M^Tot ^i^Max will 
be respectively a factor of 5 and 2 higher with respect to clusters 
which are younger and/or are allowed to form stars until recent 
times (i.e. allowing a CR with t2 = 10^ yrs). The explanation 
follows directly from the argument made for the N^^^^^ behavior 
above; matching the number of Ks -detected stars in a relatively 
old cluster with intrinsically fainter stars will require that stars 
will have to be on average more massive objects, and this will 
clearly result also in a higher total stellar mass. 

Going from IMFl to IMF3 both M^^^^ and M^Max will sig- 
nificantly increase, as expected. The trend of M^^Tot with cluster 
age is less pronounced because with IMF2 and IMF3 it is sta- 
tistically more likely to produce relatively more massive (and 
hence more easily detectable in Kg) stars requiring a lower num- 
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Fig. 11. KLF (using the absolute Ks magnitudes) for Mol3 pre- 
dicted by SCG for a CR cluster with ti = 10^ yrs and 12=10"^ yrs, 
for three diff'erent choices of the IMF (line styles as indicated 
in the figure). The dash-dotted line represents the completeness 
limit for this source given the Ks limiting magnitude. 

The shape of the resulting KLF changes throughout the Mk 
range; going from Kroupa et al.'s IMFl to Salpeter's IMF3 the 
distribution gets more skewed toward lower magnitudes; this 
was expected since IMFl produces more lower mass stars than 
IMF3. One can certainly argue that the change is not dramatic, 
but on the other hand the modification does not aff'ect one or two 
bins but the entire KLF consistently. The change is more appar- 
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ent in the region between the peak and the completeness limit 
than at the bright end of the KLF, and for this reason the ability 
of the model to discriminate among different IMFs will be better 
for those sources, as Mol3 in the figure, where the KLF's peak 
is clearly detected above the completeness limit. 
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Fig. 12. KLF (using the absolute Kg magnitudes) for Mol3 pre- 
dicted by SCG for an IMFS cluster with diff'erent ages as in- 
dicated in the plot. Older clusters produce a KLF peaked tower 
lower magnitudes. The dash-dotted line represents the complete- 
ness limit for this source given the Kg limiting magnitude. 

The diff'erence in predicted KLFs is much more dramatic if 
diff'erent age ranges are assumed, keeping fixed the shape of the 
SFH and the IMF, as it is apparent in Fig. [12] The peak of the 
KLF shifts considerably toward higher magnitudes as the me- 
dian stellar ages (t^) increase. A similar trend is observed com- 
paring SB models with diff'erent ages, although SB models al- 
ways produce KLFs which are considerably narrower than CR 
or GR models. Larger cluster ages would shift the peak of the 
KLF beyond the completeness limit; in other words, our anal- 
ysis is not sensitive to ages for the majority of stars in excess 
of ~ 5 10^ 10^ yrs; besides, such old clusters would be hard 
to justify given the fact that they are still heavily embedded in 
dense clumps. 

Finally, we briefly show how the KLFs change for difl'er- 
ent choices of the SFHs. Figure 13 shows the KLFs obtained 
for a SB with ti = 10^ yrs, compared with a CR with ti = 10^ yrs 
and t2=10^ yrs and a GR with the same start and end star for- 
mation period, and with a peak times for star formation rate of 
tc=10^ yrs. The KLFs are clearly difl'erent, with a peak magni- 
tude which is quite sensitive to the formation rate typology and 
the peak time for star production. 

As concerns the [H-K^] color functions, they are found in- 
sensitive to the choice of IMF. Similarly to the case of the KLFs, 
the main difl'erences between synthetic HKCFs are more evident 
for difl'erent age ranges especially in the number of detectable 
stars. 
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Fig. 13. KLF (using the absolute Ks magnitudes) for Mol3 pre- 
dicted by SCG for an IMFS SB cluster with ti = 10^ yrs (full fine), 
CR with ti = 10^ yrs and 12=10"^ yrs (dotted line), and GR with 
same ti and t2 and 1^=10^ yrs (dashed line). The dash-dotted line 
represents the completeness limit for this source given the Kg 
limiting magnitude. 



and where submm information was available to allow meaning- 
ful estimates of extinction (this excludes Moll 5 and Mol99). 
The number of clusters fulfilling these criteria were 16 out of 
23 detected clusters. The comparison of the observed KLFs and 
HKCFs (KLFo^^, HKCFo/?^), with the synthetic ones produced 
by SCG (KLF^^;^ HKCF^^;^), for the full set of input parameters 
(IMF, SFH and age parameters) was carried out automatically; 
to ease the process, the observed and synthetic functions were 
computed on the same M^ and H-K grid. 

The comparison procedure between synthetic and observed 
KLFs is described below, but it is the same for HKCFs. The 
KLFs are first compared bin by bin (the comparison being lim- 
ited to those bins brighter than the completeness limit) identify- 
ing with / each bin of the observed KLF, starting from i=l for 
the lower-M/^ non-zero bin to N for the bin where the complete- 
ness limit for that source is reached (the number N will clearly 
be diff'erent for each cluster). In the case of HKCFs, we exclude 
objects with H and K magnitudes greater than the observed lim- 
iting magnitudes for these bands. A matching flag m/ is set to 1 
for those bins where the number of sources coincide within the 
Icr Poissonian error bar of the observed KLF, i.e.: 



\Ni,Syn,-N^Obs] < ^JN^Obsi 



(1) 



The total number of bins where a match is found is divided 
by the total number of bins useful for the comparison to get a 
KLF compatibility figure (in %) as 



C = 100 X 



N 



(2) 



4.3. Comparing Observed and Synthetic KLFs and HKCFs 

The detailed comparison of the model KLFs and HKCFs func- 
tions to those observed was carried out only for those sources 
where the number of detected stars was sufficient (1^ >15) to al- 
low a statistically significant comparison (see Col. 3 of Table[2]), 



The higher is C, the better is the overall match between 
KLFo^^ andKLF53;„. 

However, the same value of C may result from bins con- 
centrated in the low-M^ end of the KLF, where there are few 
sources, or in the region around the peak and in the proximity of 
the completeness limit, where instead there are more sources and 
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Table 3. Results for SCG runs on detected clusters 
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hence higher statistical significance. A further figure of merit is 
then introduced, 



W . 



z 

(=1 



(3) 



Obsto 



where Ni^obstot the total number of sources present in all the 
bins used for comparison. This parameter weights each matching 
bin by its relative richness, favoring the bins closer to the com- 
pleteness limit and the KLF peak over those in the bright tail of 
the KLF, and favoring for HKCF the bin closer to the peak of the 
distribution. This choice is due to the fact that the KLF (HKCF) 
peak region is the the most sensitive to changes in the SCG input 
parameters. 

In this automatic procedure we select those models for which 
the parameters C and W (at the same time for KLFs and HKCFs) 
are maximum. For MolSB, Mol45 and Mol84 the observed KLF 
has a very irregular and multiple-peaked shape which cannot be 
matched by any model, and are therefore discarded from fur- 
ther considerations. We are then left with 14 clusters for which 
a series of models can be found with at least 75% of the bins 
matching the observations. We find that the best values of C and 
W are never found for one single set of parameters, but rather 
we identify ranges of parameter values which produce the best 
match; in other words there is a level of degeneracy that the mod- 
els cannot remove, and this varies from source to source. In 4 
clusters (Moll09, 110, 136 and 151) this degeneracy is essen- 
tially complete and the model is not able to make any predic- 
tion; in one case (Moll48) the comparison selects models with 
very old stellar ages but with total stellar luminosities by far in 
excess of the measured bolometric luminosity obtained integrat- 
ing the integrated observed for this region from the mid-IR to 
the millimeter (see Tab. In the 9 remaining case s some de- 
generacy persists especially in the IMF, confirming ( ^4.2.5 ) that 
our models are weakly selective on the IMFs, but there are clear 
indications concerning the SFH and ages. 

Table |3]reports a summary of the results. The IMF of match- 
ing classes of synthetic cluster models is shown in Col. 2. 
Cols. 3, 4 and 5 contain the times for the formation of 10%, 50% 
and 90% of the total number of cluster members; these values 
are obtained as the median of the values that these times have 
in all models that match the observations. Col. 6 is the number 
of cluster members N^; Col. 7 shows the mass of the most mas- 
sive object M^Max'^ Cols. 8 reports the total stellar mass of the 
cluster Mi,Tot- We stress again that the analysis selects classes of 
models rather than single models; the values reported in Table |3] 
are the median of the values that the parameters assume within 
the class of matching models. The table shows that for some 



fields multiple IMFs are compatible with the data and, in gen- 
eral, SFHs with constant (CR) or Gaussian (GR) star formation 
rate provide acceptable solutions for certain age ranges (as re- 
ported in the Table). Simulations with SFHs with a single burst 
are, in general, not accepted. Our modeling is insensitiv e to bu lk 
stellar ages in a cluster in excess of 5 10^ ^ 10^ years (^4.2.5). 



5. Discussion 

5.1. Cluster Ages and Star Formation Histories 

Perhaps the most important result of this work is that in all clus- 
ters where the comparison of observed KLFs with the ones pre- 
dicted by the SCG model is possible (see previous paragraph), 
the observations are consistent with a star formation which goes 
on over time intervals that in most cases have a spread between 
about few 10^ and few 10^ yrs , and with a median cluster age 
of a few 10^ yrs. In most cases we cannot discriminate clearly 
between a constant or variable SFR but we are confident that 
we can exclude that on average the stars in our clusters are co- 
eval and originating from a single burst of formation. Detailed 
studies toward the Orion Nebula Cluster show that stars have 
been forming for at least 1 tdyn, or 20 i f f (Pall a & Stabler 19991 
Hoogerwerf et al. 120011 Hillenbrand I1997K and our results 
would seem to generalize this on a larger sample of intermediate 
and high-mass star forming regions. 

In principle it can be argued that our analysis is incomplete 
since we did not take in consideration longer wavelength data 
which could pinpoint heavily extincted objects barely visible, or 
not visible at all, in the near infrared. This however, does not 
appreciably modify our conclusions about the age spread within 
the clusters. Indeed, Vig et al. (12007 ) applied a diff'erent analy- 
sis to the specific region Mol075, a field not included in our fi- 
nal analysis (Table|3]) because the background- subtracted KLF is 
populated by too few objects for a statistically significant model 
comparison. Vig et al. also considered Spitzer IRAC and MIPS 
data, looking for the brighter and redder objects in the area cov- 
ered by submillimeter emission. In this way they could identify 
the younger and more massive objects in the field with an es- 
timated age of the order of 10^ years or less. This approach, 
however, is not sensitive to low mass and relatively older pre- 
MS objects, for which our method is ideally designed. While for 
this particular field, for the reasons explained above, we cannot 
perform a direct comparison to our approach, it is clear that the 
inclusion of longer wavelength data in the analysis might have 
identified a diff'erent, younger, population of objects, rather am- 
plifying the observed age spread deduced for the clusters. 

Models of cluster formation via competitive accretion seem- 
ingly succeed in delivering an IMF close to the observed ones 
thanks to the spread of the accretion rates consequent to the com- 
petitive accretion mechanism, but the prediction that all stars are 
formed in about 5x10^ yrs (Bonnel et al. 2004 ) for typical condi- 
tions in young clusters, corresponding to a dynamical time or so, 
seems to be in disagreement with our results. We instead favour 
scenarios (Tan & McKee 170021) where stars keep forming over 
several free-fall times thus providing the required age spread. 
The finding that the most massive object in the fields considered 
in this work are still being formed or have just finished a phase 
of intense accretion (Molinari et al. 2008a) is a further indica- 
tion that star formation seems to be a long-duration process in 
the life of a molecular clump. 

How do our clusters compare to relatively more evolved sys- 
tem, like those observed for a sample of Herbig Ae/Be by Testi et 
al. (.1997til998i) ? Figure [14] shows the relationship between the 
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Fig. 14. Number of cluster members as a function of Mass of highest mass star. Asterisks are for Testi et al.'s (|1997|[T998l ) Herbig 
Ae/Be sample. Other point are our source sample, where the full diamonds are the obtained from our clusters simulation analysis; 
the full lines passing through the full diamonds represent the total spread of the parameters from all SCG models that match the 
observed KLF and HKCF. Full triangles are the same clusters where this time we use the observationally derived 1^ instead of 
the model-provided N^. Empty triangles are those clusters which were not analyzed with SCG, or which exhibit complete model 
degeneracy; in these cases Ic was used for the cluster membership (from Table [2]), while the maximum stellar mass has been derived 
assuming that half of the bolometric luminosity is generated by a single ZAMS star (the horizontal lines through the empty triangles 
represent the mass spread assuming that a fraction between 30% and 100% of L^^/ come from a single star) 



mass of the most massive source and the total number of cluster 
stars as provided by the SCG simulations for our modeled clus- 
ters. The filled symbols represent the clusters which we could 
model (Tab.[3|; the empty symbols represent the clusters which 
instead could not be modeled for a variety of reasons (see ^4.3[ ), 
while Testi et al.'s clusters are reported as asterisks (see figure 
caption for detailed explanation of the symbols). The figure sug- 
gests that the clusters presented in this study are richer than those 
surrounding Herbig Ae/Be stars for any given value of the most 
massive star in each cluster. The trend persists if we use similar 
indicators (e.g. 1^, the full triangles in the figure). Furthermore, 
we note that while the limiting magnitudes of our observations 
and those of Testi et al. ( 1998 ) are similar, higher Ay values to- 
ward our sources and the typically greater distance from the Sun 
would justify the non detection of the fainter cluster members 
predicted by the SCG models. It is thus likely that the values of 
Ic derived from our observations tend to underestimate the clus- 
ter memberships. 



This evidence is clear for values of the most massive star in 
the cluster below ~ lOM©, where the 9 clusters for which we 
could compare observations with SCG predictions lie (the dia- 
monds). In the clusters for which this could not be done (the 
empty triangles), an estimate of the mass for the highest-mass 
star was made by assuming that a fraction between 30 and 100% 
of the bolometric luminosity is due to a single ZAMS star. In this 
case the trend for richer clusters compared to Herbig stars (i.e., 
the asterisks) would tend to become marginal. These estimates, 
however, place the latter clusters systematically to the right in 
the plot, compared to the 9 modeled clusters; indeed, if we were 
to estimate in the same way a maximum stellar mass also for 
the 9 modeled clusters, we would obtain values in excess (be- 
tween a factor two and three) to those provided by the detailed 
SCG modeling. In other words, the evidence that our clusters are 
richer than those around Herbig stars is marginal at worst (i.e., 
using the most conservative approach of estimating the mass of 
the highest-mass star). 
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Fig. 15. Distribution of radii for our clusters (full line) and those 
associated with Herbig Ae/Be stars (from Testi et al. (2001 K 
dashed line). 



The plausibility of this interpretation is strengthened by the 
recent results of Baumgardt & Kroupa ( 2007 ) who carried out 
extensive numerical simulations of the evolution of stellar clus- 
ters as a function of, among other parameters, the cluster gas 
content. They show that for a wide range of initial conditions 
and star formation efficiencies, the dispersal of the gas with age 
causes the cluster to expand overall and disperse a fraction of 
the stars originally belonging to the cluster. Besides, as the clus- 
ter expands its decreasing stellar density would make the low- 
density outer regions of the cluster ever more difficult to detect 
against the field stars (especially in the Galactic plane, where all 
these objects lie), mimicking a smaller cluster from an observa- 
tional viewpoint. Fig. [15] shows the distribution of radii for our 



clusters (full line) and for those surrounding Ae/Be stars; radii 
have been derived with the same analysis in the two samples and 
the histogram clearly shows that our clusters are indeed larger in 
size, confirming the prediction of Baumgardt & Kroupa ( 2007 ). 

This age eff'ect on clusters size is also revealed within our 
clusters sample. Fig. 16 presents the relationship between the 
clusters radii as derived from the observations and their ages 
as derived from the modeling. The reported correlation has a 
Spearman rank correlation coefficient of ~ -0.6 indicating a 
good correlation, with a significance of 92% (between 2 and 
3cr). Ongoing gas dispersal is certainly plausible in our clus- 
ters, given the common detection in these systems of molecular 
outflows (Zhang et al. 12001 1 [2005 b which are highly efl'ective in 
transferring material away from the forming objects and possi- 
bly out of the star forming region; parsec-scale flows are indeed 
commonly observed also from low-mass YSOs. 

The final stage of gas dispersal, eventually leading to opti- 
cally revealed clusters like those around Ae/Be stars, might be 
triggered by the powerful winds and radiation fields from new- 
born O and B stars; indications are (Molinari et al. 2008a, 2008b ) 
that the most massive objects forming in the densest regions of 
the clumps hosting our clusters may have not yet reached the 
ZAMS, or are just starting to develop their Hii regions. It is quite 
likely that this event will mark the moment of maximum effi- 
ciency of gas dispersal and further evolution of our clusters to- 
ward the Ae/Be 's ones. 




Fig. 16. Cluster radii (from Tab. [2| as a function of the cluster 
median ages (from Tab. [3]). Dashed lines represent the linear fits 
obtained fitting in turn one plotted variable as a function of the 
other; the full line is the bisector of the two dashed lines and 
represents the fit which minimizes the quadratic geometric (i.e. 
not along the X or Y axis alone) distance of the data from the fit. 
Following Isobe et al. (1990) , this is the adequate approach when 
the nature of the data scatter around the linear fit is not known 
(and it is not due to classical measurement uncertainties); the 
slope of the full line is -0.8 ± 0.2 and the Icr spread is within the 
area enclosed by the two dashed lines. 



5.2. Physical vs. Statistical Models for Cluster Formation 

Figure [14] can also be used as a diagnostic to discriminate be- 
tween diff'erent classes of models for the origin of clusters. Testi 
et al. (1200 ll) called physical the class of models which imply a 
physical relationship between the most massive star that forms 
and other environmental properties like the cluster richness or 
the mass of the gaseous clump where the stars originate from; 
examples are the "turbulent core" (McKee & Tan 2003 ), the "co- 
alescence" (Bonnell et al., 1998), or the competitive accretion 
models (Bonnell & Bate, 2006). In a second class of models, 
called statistical, the relationship between the most massive star 
in a cluster and its richness arises from the higher probability 
of finding the rare massive stars in rich clusters rather than in 
isolation (Bonnell & Clarke 1 19991) . If clusters are populated by 
randomly picking stars from the field stars' initial mass func- 
tion, and considering a cluster membership-size distribution in 
the form of an appropriate power law, then the observations of 
Testi et al. (1999 ) can be naturally explained. Nevertheless, this 
model predicts that a significant fraction of high-mass stars are 
still associated with relatively poorly populated clusters, in other 
words that massive stars can be found both in high-N clusters 
and, to a lesser extent, in low-N clusters. 
The dashed line in Fig. 



14 is the upper boundary of the re- 



gion which should contain 25% of the cluster realizations ob- 
tained by randomly extracting stars from the IMF (Testi et al. 
2001 ). Indeed, if we consider our measurements of 1^ for our en- 
tire sample (i.e. the full and empty triangles), there is fraction of 
about 15% of the clusters which is found marginally below the 
dashed line. However, it has to be noted that our modeling was 
possible only for clusters above a membership threshold derived 
with Ic, it is thus a somewhat biased subsample toward rich clus- 
ters. In the extreme assumption that the fields with no detected 
cluster are cases of systems composed by a single heavily ex- 
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tincted star, and thus would fall below the dashed line in Fig. 14 
then this fraction would approach the 25% level. This, however, 
is an extreme case because, as we already discussed, the high 
value of the extinction derived from submillimeter maps may be 
the reason for not detecting clusters in at least a number of ob 
served fields. 



> 1.0 




Fig. 17. Mass for the highest mass star as a function of the total 
stellar mass for the 9 modeled clusters (see Table |3j; the bars 
associated to each cluster (the diamond symbols) represent the 
total span of the parameter values for the classes of models se- 
lected by our analysis (the diamond marks the median values). 
The dash-dotted, and dashed lines represents the relationship ob- 
tained for numerical simulations of clusters drawn from pure 
random sampling of the IMF, and using a so-called sorted sam- 
pling, following Weidner & Kroupa (2006). The full line is a 
semi-analytical approximation of this relationship obtained by 
Weidner & Kroupa ( 2004 ). The dotted line is the limit where a 
cluster is made of just one star. 



As a further diagnostic between physical and statistical clus- 
ter models, Weidner & Kroupa (2006) recently argued that a 
non-trivial correlation exists between the highest-mass star in 
a cluster, M^uax^ and the total stellar mass of the cluster itself. 



(Fig. 17). Numerical simulations show that the relation- 
ship obtained by pure random sampling of the IMF with an im- 
posed physical limit of 150 M© for the maximum stellar mass 
(the dashed-dot line in Fig. 17 ) is clearly not representing our 
data. 

A substantially diff'erent result (the dashed line) is obtained 
if cluster members are picked in ascending order and constrained 
to total cluster masses distributed according to a cluster total 
mass function (Weidner & Kroupa 2006). Basically, this second 
option represents the fact that drawing 10 cluster of 100 stars 
will not deliver the same M^Max^^Q^i^Tot) as drawing 1 cluster 
of 1000 stars. This trend closely resembles a semi-analytical ap- 
proximation of M^Mflx=f(M*ro?) obtained by Weidner & Kroupa 
((2004), again assuming that total cluster stellar masses are dis- 
tributed according to a power-law mass function. Weidner & 
Kroupa (2006 ) suggest that this sorted sampling way of pop- 
ulating a cluster can be physically understood in terms of a 
pre-stellar clump where initial low-amplitude perturbations start 
low-mass star formation; as further perturbations with larger am- 
plitude grow, higher mass stars will start to form until the feed- 
back from the latter will start disrupting the cloud. This scenario 



of an organised star formation where low-mass stars are the first 
ones to form, is the same we favor (see %5.\\ considering the 
age spread we find in our clusters in which, based on indepen- 
dent considerations (Molinari et al. 2008a), the most massive star 
may have not been formed at all. By the way, this latter possi- 
bility does not change the substance of the agreement between 
our data and the physical cluster models in Fig. 17 since the late 
addition of the highest mass star currently not yet visible in the 
near-IR would shift the points toward the top-right of the plot. 

We verified a posteriori that the range of M^m^x and M^^Tot 
parameters values explored by our simulations is much wider 
than the area spanned by the bars attached to the single points, 
and also includes regions compatible with the random sampling 
cluster model. We then conclude that the result of Fig.fTTlis not 



likely to be produced by a biased sampling of the clusters' phys- 
ical parameters explored in our models. 

The predictions of the sorted sampling above descibed, with 
which our data points best agree, are also in good agreement 
with the results from simulations of clusters forming in a com- 
petitive accretion scenario (Bonnell, Vine & Bate 120041) . This 
model, however, seems excluded by the observed ages and age 
spreads in our clusters which are in clear disagreement with the 
predicted cluster formation timescales of 2-^3 free-fall times. 

5.3. Influence of Binarity on the Interpretation of Age Spread 

Weidner, Kroupa & Maschberger (2008 ) carried out extensive 
numerical simulations to determine how the presence of unre- 
solved binary/multiple stars can aff'ect the observational prop- 
erties of a young cluster in a massive star forming region. 
Assuming a 100% of binarity in a cluster of coeval sources, they 
find that unresolved binaries may lead an observer to conclude 
that instead a significant age spread is present in the cluster; the 
full line in Fig. [18] shows the a posteriori age determination as- 
suming that all binaries are unresolved. We see that the measured 
age spread for the large majority of stars simulated in Weidner 
et al.'s simulation is comparable to a Log(age) Gaussian distri- 
bution with (T=0.1, which is one of the possible choices of Star 
Formation Histories in our cluster models. However, the com- 
parison of our observed KLFs with the SCG models (Table [3]) 
suggests age spreads much bigger than this, and more compara- 
ble to a Log(age) distribution with (T=0.5 (the dotted line in Fig. 
[18]). We then conclude that unresolved binaries in our clusters 
cannot account for the observed age spread. 

6. Conclusions 

The main results of this work are the following: 

- We have imaged in the J, H, and Ks NIR bands 26 inter- 
mediate and high-mass star forming regions selected from a 
larger sample of sources and spanning a range in luminosi- 
ties and presumed youth. We have identified the presence of 
23 young stellar clusters in 22 fields. 

- Revealed clusters have richness indicator values between ten 
and several tens of objects and have median radius of 0.7 pc. 
Compared to clusters around Herbig Ae/Be stars, our clus- 
ters seem richer and larger for any given mass for the most 
massive star in each cluster. Color-color diagrams show that 
these clusters are young: many sources show colors typical 
of young pre-MS objects with an intrinsic IR excess arising 
from warm circumstellar dust. This is confirmed by the anal- 
ysis of color-magnitude diagrams where a significant frac- 
tion of stars in each cluster are found in correspondence to 
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Fig. 18. Full line represents the age spread resulting from the 
simulations of Weidner, Kroupa and Maschberger (2008 ) for a 
cluster with an input age of 2 x 10^ yrs and 100% binarity frac- 
tion. The two gaussians are the age weight functions used in our 
simulations of Gaussian Star Formation Histories, with cr=0.l 
(dashed line) and 0.5 (dotted line) respectively. In case of con- 
stant SFH models, the adopted age spread is comparable to the 
cr=0.5 distribution in the figure. 



the Pre-MS evolutionary tracks even after conservative de- 
reddening is applied. 

- We cannot perform a direct inversion of stellar luminosities 
(and colors) into masses and ages; we use a Synthetic clus- 
ter Generator (SCG) model to create statistically- significant 
cluster simulations for diff'erent initial parameters (IMF, 
SFH, source ages and their distribution), and compare the 
synthetic KLFs and HKCFs with the observed (field star- 
subtracted) ones. For the fraction of clusters for which this 
comparison selects a well-confined region of the parameter 
space, we conclude that star formation in these regions can- 
not be represented with a single burst, but is a process that is 
spread out in time. Clusters have mean ages of a few 10^ yrs; 
the ages of most of the clusters members are spread, within 
each cluster, between a few 10^ yrs to a few 10^ yrs. Together 
with the independent evidence that the most massive stars in 
these systems are very young, or not even on the ZAMS yet, 
this result is difficult to reconcile with any model predicting 
cluster formation in a crossing time. 

- The cluster radii seem to be inversely proportional to their 
age, as also confirmed by the comparison of cluster parame- 
ters with those typical of Ae/Be systems, which are smaller 
and less rich. As suggested by numerical simulations in the 
literature, dispersal of intra-cluster gas (by, e.g., molecular 
outflows or radiation fields from massive stars) may lead to 
the loss of a fraction of cluster stellar population, thus indeed 
leading to smaller and less rich clusters. Our results seem in 
line with this prediction. 

- The relation between the mass of the most massive star in 
a cluster and the cluster's richness indicator suggests that a 
physical rather than statistical nature of the cluster origin is 
more likely. 
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Appendix A: images for all observed fields 

This appendix contains the images for all ob served 
fields; these are also available at http://galate a.ifsi- 
roma.inaf.it/faustini/K-images/ 
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00''45"10' 00' 60' 

a (J2000) 



Fig. A.l. image of field Mol003 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 1200 8 at . 




a (J2000) 



Fig. A.2. image of field MolOOS with superimposed SCUBA 850yum continuum in white contours (Molinari et al. l2008al) . 
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Fig. A.3. image of field Mol009; no submillimeter image is available for this field. 
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Fig. A.4. image of field MolOl 1 with superimposed SCUBA 850yL/m continuum in white contours (Molinari et al. l2008al) . 
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5''40"30' 25* 20* 

a (J2000) 



Fig. A.5. image of field Mol012 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 1200 8 at . 




6''08~40' 
a (J2000) 



Fig. A.6. image of field Mol015; no submillimeter image is available for this field. 
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a (J2000) 



Fig. A.7. image of field Mol028; no submillimeter image is available for this field. 




a (J2000) 



Fig. A.8. image of field Mol030; no submillimeter image is available for this field. The white cross marks the position of the 
IRAS source. 
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18'*05"30' 25* 20* 

a (J2000) 



Fig. A.9. image of field Mol038 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 1200 8 at . 




a (J2000) 



Fig. A.IO. image of field Mol045 with superimposed SCUBA 850//m continuum in white contours (Molinari et al. 12008 at . 
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18''19"10' 05* 
a (J2000) 



Fig. A.ll. image of field Mol050; no submillimeter image is available for this field. The white cross marks the position of the 
IRAS source. 




a (J2000) 



Fig. A.12. image of field Mol059 with superimposed SCUBA 850//m continuum in white contours (Molinari et al. l2QQ8ab . 
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a (J2000) 



Fig. A.13. image of field Mol075 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 12008 at . 




a (J2000) 



Fig. A.14. image of field Mol082; no submillimeter image is available for this field. The white cross marks the position of the 
IRAS source. 
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a (J2000) 



Fig. A.15. image of field Mol084; no submillimeter image is available for this field. The white cross marks the position of the 
IRAS source. 




19''ll"45' 40* 35' 

a (J2000) 



Fig. A.16. image of field Mol098 with superimposed SCUBA 850//m continuum in white contours (Molinari et al. l2QQ8ab . 
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a (J2000) 



Fig. A.17. image of field Mol099; no submillimeter image is available for this field. The white cross marks the position of the 
IRAS source. 




a (J2000) 



Fig. A.18. image of field Moll03 with superimposed SIMBA 1.2mm continuum in white contours ( Beltran et al. l2QQ6l) . 
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19"'39"36* 30* 
a (J2000) 



Fig. A.19. image of field Moll09 with superimposed SIMBA 1.2mm continuum in white contours ( Beltran et al. l2006t . 




a (J2000) 



Fig. A.20. image of field Moll 10 with superimposed SIMBA 1.2mm continuum in white contours ( Beltran et al. l2006t . 
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2l''32"40' 30' 
a (J2000) 



Fig. A.21. image of field Moll 36 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 12008 at . 




a (J2000) 



Fig. A.22. image of field Moll 39 with superimposed SCUBA 850//m continuum in white contours (Molinari et al. 12008 at . 
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Fig. A.24. image of field Moll48 with superimposed MAMBO 1.1mm continuum in white contours (Molinari et al. l2008al) . 
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a (J2000) 



Fig. A.25. image of field Moll 51 with superimposed SCUBA 850yum continuum in white contours (Molinari et al. 12008 at . 




a (J2000) 



Fig. A.26. image of field Moll 60 with superimposed SCUBA 850//m continuum in white contours (Molinari et al. 12008 at . 
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Appendix B: luminosity functions 

This appendix presents the set of background- subtracted lu- 
minosity functions for all detected clusters; the green vertical 
line represents the 90% completeness limit as estimated from 
artificial star experiements (see ^2.1|). Material can also be re- 
trieved at http://g alatea.ifsi-roma.inaf.it/faustini/KLF/^ 
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Fig.B.l. M0IOO3 (left), M0IOO8A (center), M0IOO8B (right) 




Absolute Magnitudes Absolute Magnitudes 

Fig.B.2. M0IOO9 (left), MolOll (center), Mol012 (right) 
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Fig.B.3. M0IOI5 (left), M0IO28 (center), Mol045 (right) 
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Absolute Magnitudes 



Fig.B.4. M0IO5O (left), M0IO75 (center), Mol082 (right) 
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Fig.B.5. Mol084 (left), Mol099 (center), Moll03 (right) 
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Fig.B.6. Moll36 (left), Moll 10 (center), Moll36 (right) 
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Fig.B.7. M0II39 (left), M0II43 (center), Moll48 (right) 
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Fig.B.8. M0II5I (left), M0II6O (right) 
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